library(coda)
library(bayesplot)
library(ggplot2)
library(ggsci)
library(khroma)
library(tidyverse)
library(reshape2)
library(yardstick)
library(here)
knitr::opts_chunk$set(echo = TRUE, dpi = 300 )
Set up MCSim file
# this markdown file must be saved in top level directory for the following to work; the mcsim code depends on getwd results.
mdir <- "MCSim"
source(here::here(mdir,"setup_MCSim.R"))
# Make mod.exe (used to create mcsim executable from model file)
makemod()
The mod.exe had been created.
id_lut <- multicheck$df_check %>% select(Level) %>% unique () %>%
mutate(dataset = c(
rep("Decatur M Train", 9),
rep("Decatur F Train", 9),
rep("Decatur M Test", 9),
rep("Decatur F Test", 10),
'Paulsboro-Train','Horsham-Train',
'Warminster-Test'),
Sex = c(
rep("M", 9),
rep("F", 9),
rep("M", 9),
rep("F", 10),
'Mixed', 'Mixed', 'Mixed'),
City = c(
rep("Decatur", 18),
rep("Decatur", 19),
'Paulsboro','Horsham','Warminster'),
Train_Test = c(
rep("Train", 9),
rep("Train", 9),
rep("Test", 9),
rep("Test", 10),
'Train', 'Train', 'Test'),
datatype = c(
rep("Individual",9+9+9+10),
rep("Summary",3)),
Simulation = row_number(),
variable = paste0(dataset, " ",Simulation))
id_lut$dataset <- factor(id_lut$dataset,levels=
c("Decatur M Train","Decatur F Train","Arnsberg M Train",
"Arnsberg F Train","Decatur M Test","Decatur F Test","Arnsberg M Test",
"Arnsberg F Test","Minnesota Train","Minnesota Test",
'Lubeck-Bartell-Train', 'Lubeck-Bartell-Test',
'Little Hocking-Bartell-Train', 'Little Hocking-Bartell-Test',
'Little Hocking-Emmett-Test','Paulsboro-Train','Horsham-Train',
'Warminster-Test','Warrington-Train'))
id_lut$City <- factor(id_lut$City,levels =
c("Decatur","Arnsberg","Minnesota",'Lubeck-Bartell',
'Little Hocking-Bartell','Little Hocking-Emmett',
'Paulsboro','Horsham','Warminster','Warrington'))
indiv_lut <- id_lut %>%
filter(City %in% c("Decatur")) %>%
mutate( dataset = as.factor(dataset))
nv <- data.frame(dataset =unique(indiv_lut$dataset),
variable= rep("Pop GM", 4),
type= rep("Pop GM", 4), stringsAsFactors = FALSE)
set.seed(314159)
indiv_parms <- indiv_lut
lnkparmnames <- paste("ln_k.",gsub("_",".",indiv_parms$Level),".",sep="")
lnVdparmnames <- paste("ln_Vd.",gsub("_",".",indiv_parms$Level),".",sep="")
parmsamp <- apply(multicheck$parms.samp,2,sample,1)
## Random z-score estimate of each parameter
indiv_parms$ln_k.z.samp <- parmsamp[lnkparmnames]
indiv_parms$ln_Vd.z.samp <- parmsamp[lnVdparmnames]
normd <- data.frame(x=qnorm(ppoints(200)))
normd$y <- dnorm(normd$x)
iplotk<-
ggplot(subset(indiv_parms,Train_Test=="Train"))+
geom_histogram(aes(x=ln_k.z.samp,after_stat(density)),bins=10)+facet_wrap(~City,ncol=1)+
geom_line(aes(x=x,y=y),data=normd)+
xlab("Individual z-scores for k") + theme_bw()
iplotVd<-
ggplot(subset(indiv_parms,Train_Test=="Train"))+
geom_histogram(aes(x=ln_Vd.z.samp,after_stat(density)),bins=10)+facet_wrap(~City,ncol=1)+
geom_line(aes(x=x,y=y),data=normd)+
xlab("Individual z-scores for Vd") + theme_bw()
print(iplotk)
print(iplotVd)
ggsave(file.path("output-plots",
paste0( sa,"Indiv_zscores_k_PFNA.pdf")),iplotk,dpi=600)
Saving 3.5 x 3.5 in image
ggsave(file.path("output-plots",
paste0( sa,"Indiv_zscores_Vd_PFNA.pdf")),iplotVd,dpi=600)
Saving 3.5 x 3.5 in image
ggsave(file.path("output-plots",
paste0( sa,"Indiv_zscores_k_PFNA.png")),iplotk,dpi=600)
Saving 3.5 x 3.5 in image
ggsave(file.path("output-plots",
paste0( sa,"Indiv_zscores_Vd_PFNA.png")),iplotVd,dpi=600)
Saving 3.5 x 3.5 in image
This is a Figure 2 panel. Needed to use “scale=1.1” in ggsave to match PFOA.
nrow(multicheck$df_check)
[1] 38500
nrow(id_lut)
[1] 40
multicheck$df_check %>% left_join(id_lut) %>% nrow()
Joining, by = c("Level", "Simulation")
[1] 38500
names(multicheck$df_check)
[1] "Level" "Simulation" "Output_Var" "Time" "Data"
[6] "Prediction"
Level
Simulation
Output_Var
Time
Data
Prediction
multicheck2 <- multicheck$df_check %>% left_join(id_lut)%>%
group_by_at ( vars(-Prediction)) %>%
summarise(Prediction = median(Prediction)) %>%
ungroup() %>%
group_by(City) %>%
mutate(Train_Test = factor(Train_Test, levels = c("Train", "Test")),
`City (datatype)` = factor (paste0("\n" , City, "\n(", datatype, ")\n") ),
label = case_when(Train_Test=="Train" ~ "E: PFNA Train",
Train_Test=="Test" ~ "F: PFNA Test",
TRUE ~ ""))
Joining, by = c("Level", "Simulation")
`summarise()` has grouped output by 'Level', 'Simulation', 'Output_Var', 'Time', 'Data', 'dataset', 'Sex', 'City', 'Train_Test', 'datatype'. You can override using the `.groups` argument.
#define color for testing boxplots
bp_cols <- c (as.character (khroma::colour("muted")(9)) , "#191919")
bp_cols <-bp_cols[c(1, 7, 9:8)]# plot_scheme_colourblind(bp_cols)
### Create aesthetics lookup
aes_lut <- multicheck2 %>% ungroup() %>%
group_by(City, datatype, `City (datatype)` ) %>% summarise () %>% ungroup() %>%
mutate( cols = bp_cols, city_fills = bp_cols ,
# for individual level on point plot (multicheck2), darken outlines for visibility, use standard colors otherwise
city_outlines = if_else(datatype == "Individual" , colorspace::darken(city_fills, 0.3), city_fills) ,
shapes = case_when(datatype == "Individual" & `City` %in% c('Decatur', 'Arnsberg', 'Minnesota') ~ 23,
datatype == "Summary" &`City` %in% c("Horsham", "Warminster", "Warrington") ~ 2,
datatype == "Summary" & `City` == "Paulsboro" ~ 1,
TRUE ~ 18 ),
size = if_else(datatype =="Individual", 1.75, 2.5 ) )
`summarise()` has grouped output by 'City', 'datatype'. You can override using the `.groups` argument.
source( paste0(gsub(basename(here()), 'shared_functions', here()), '/plot_scatter_mcheck.r'))
p2 <- plot_scatter_mcheck(dframe = multicheck2, pfas_nom = pfas_name, aes_lut_fn = aes_lut )
Joining, by = "facet_label"
print(p2)
ggsave(here ("output-plots", paste0( sa,"multicheckplot_", pfas_name,
".pdf")),p2,dpi=600, scale=1.1)
Saving 8.8 x 3.85 in image
ggsave(here ("output-plots", paste0( sa,"multicheckplot_", pfas_name,
".png")),p2,dpi=600, scale=1.1) # doesn't properly handle '\n' unicode for line feed. Not sure why; possibly device opens from windows working environment.
Saving 8.8 x 3.85 in image
multicheck2$PFAS <- "PFNA"
multiplot_scatter_dat <- list(aes_lut = aes_lut, multicheck2 = multicheck2) # save out for use in commbined plot
save(multiplot_scatter_dat, file= "pfna_multiplot_scatter.Rdata")
df_check <- multicheck$df_check
df_check <- subset(df_check,Data > 0)
n1 <- nrow(df_check)
id_chks <- df_check %>% select(Level) %>% unique() %>% bind_cols(id_lut) %>%
mutate(dataset = as.factor(dataset), Sex = as.factor(Sex), City = as.factor(City),
Train_Test = as.factor(Train_Test))
New names:
* Level -> Level...1
* Level -> Level...2
df_check <- df_check %>% left_join(id_chks, by = "Simulation")%>%
mutate(Dataset = paste(as.character(dataset), Simulation),
Sex = ordered(Sex, levels = c("M", "F", "Mixed"),
labels = c("Male", "Female", "Mixed (all sexes)")))
n2 <- nrow(df_check)
if(n1 != n2)print("duplicates created in id-lut join")
df_check$Time.desc <- as.character(paste0("T=",df_check$Time))
df_check$Time.desc[df_check$Time.desc == "T=1e-06"] <- "SteadyState"
df_check$Dataset.Time <- interaction(df_check$Dataset,
df_check$Time.desc,lex.order=TRUE)
df_check$Dataset.Time <- factor(df_check$Dataset.Time,
levels=levels(df_check$Dataset.Time))
calibdata <- df_check[,names(df_check) != "Prediction"]
calibdata <- calibdata[!duplicated(calibdata),]
print(calibdata)
Level Simulation Output_Var Time Data Level...1 Level...2
1 1_1_1_1_1 1 Cserum 0.000 1.900 1_1_1_1_1 1_1_1_1_1
2 1_1_1_1_1 1 Cserum 5.802 1.100 1_1_1_1_1 1_1_1_1_1
3 1_1_1_1_2 2 Cserum 0.000 1.100 1_1_1_1_2 1_1_1_1_2
4 1_1_1_1_2 2 Cserum 5.802 0.900 1_1_1_1_2 1_1_1_1_2
5 1_1_1_1_3 3 Cserum 0.000 2.100 1_1_1_1_3 1_1_1_1_3
6 1_1_1_1_3 3 Cserum 5.802 0.700 1_1_1_1_3 1_1_1_1_3
7 1_1_1_1_4 4 Cserum 0.000 1.900 1_1_1_1_4 1_1_1_1_4
8 1_1_1_1_4 4 Cserum 5.802 0.500 1_1_1_1_4 1_1_1_1_4
9 1_1_1_1_5 5 Cserum 0.000 3.500 1_1_1_1_5 1_1_1_1_5
10 1_1_1_1_5 5 Cserum 5.802 0.900 1_1_1_1_5 1_1_1_1_5
11 1_1_1_1_6 6 Cserum 0.000 1.200 1_1_1_1_6 1_1_1_1_6
12 1_1_1_1_6 6 Cserum 5.802 0.600 1_1_1_1_6 1_1_1_1_6
13 1_1_1_1_7 7 Cserum 0.000 2.700 1_1_1_1_7 1_1_1_1_7
14 1_1_1_1_7 7 Cserum 5.802 0.800 1_1_1_1_7 1_1_1_1_7
15 1_1_1_1_8 8 Cserum 0.000 1.700 1_1_1_1_8 1_1_1_1_8
16 1_1_1_1_8 8 Cserum 5.802 0.500 1_1_1_1_8 1_1_1_1_8
17 1_1_1_1_9 9 Cserum 0.000 1.100 1_1_1_1_9 1_1_1_1_9
18 1_1_1_1_9 9 Cserum 5.802 0.600 1_1_1_1_9 1_1_1_1_9
19 1_1_1_1_10 10 Cserum 0.000 1.600 1_1_1_1_10 1_1_1_1_10
20 1_1_1_1_10 10 Cserum 5.802 1.000 1_1_1_1_10 1_1_1_1_10
21 1_1_1_1_11 11 Cserum 0.000 2.000 1_1_1_1_11 1_1_1_1_11
22 1_1_1_1_11 11 Cserum 5.802 1.100 1_1_1_1_11 1_1_1_1_11
23 1_1_1_1_12 12 Cserum 0.000 1.000 1_1_1_1_12 1_1_1_1_12
24 1_1_1_1_12 12 Cserum 5.802 0.400 1_1_1_1_12 1_1_1_1_12
25 1_1_1_1_13 13 Cserum 0.000 1.700 1_1_1_1_13 1_1_1_1_13
26 1_1_1_1_13 13 Cserum 5.802 0.400 1_1_1_1_13 1_1_1_1_13
27 1_1_1_1_14 14 Cserum 0.000 3.900 1_1_1_1_14 1_1_1_1_14
28 1_1_1_1_14 14 Cserum 5.802 1.200 1_1_1_1_14 1_1_1_1_14
29 1_1_1_1_15 15 Cserum 0.000 2.400 1_1_1_1_15 1_1_1_1_15
30 1_1_1_1_15 15 Cserum 5.802 3.600 1_1_1_1_15 1_1_1_1_15
31 1_1_1_1_16 16 Cserum 0.000 1.200 1_1_1_1_16 1_1_1_1_16
32 1_1_1_1_16 16 Cserum 5.802 0.600 1_1_1_1_16 1_1_1_1_16
33 1_1_1_1_17 17 Cserum 0.000 2.900 1_1_1_1_17 1_1_1_1_17
34 1_1_1_1_17 17 Cserum 5.802 0.900 1_1_1_1_17 1_1_1_1_17
35 1_1_1_1_18 18 Cserum 0.000 1.200 1_1_1_1_18 1_1_1_1_18
36 1_1_1_1_18 18 Cserum 5.802 0.400 1_1_1_1_18 1_1_1_1_18
37 1_1_1_2_1 19 Cserum 0.000 1.500 1_1_1_2_1 1_1_1_2_1
38 1_1_1_2_1 19 Cserum 5.802 0.400 1_1_1_2_1 1_1_1_2_1
39 1_1_1_2_2 20 Cserum 0.000 3.300 1_1_1_2_2 1_1_1_2_2
40 1_1_1_2_2 20 Cserum 5.802 1.000 1_1_1_2_2 1_1_1_2_2
41 1_1_1_2_3 21 Cserum 0.000 1.800 1_1_1_2_3 1_1_1_2_3
42 1_1_1_2_3 21 Cserum 5.802 0.600 1_1_1_2_3 1_1_1_2_3
43 1_1_1_2_4 22 Cserum 0.000 1.600 1_1_1_2_4 1_1_1_2_4
44 1_1_1_2_4 22 Cserum 5.802 0.700 1_1_1_2_4 1_1_1_2_4
45 1_1_1_2_5 23 Cserum 0.000 4.300 1_1_1_2_5 1_1_1_2_5
46 1_1_1_2_5 23 Cserum 5.802 1.400 1_1_1_2_5 1_1_1_2_5
47 1_1_1_2_6 24 Cserum 0.000 1.700 1_1_1_2_6 1_1_1_2_6
48 1_1_1_2_6 24 Cserum 5.802 1.400 1_1_1_2_6 1_1_1_2_6
49 1_1_1_2_7 25 Cserum 0.000 1.300 1_1_1_2_7 1_1_1_2_7
50 1_1_1_2_7 25 Cserum 5.802 0.400 1_1_1_2_7 1_1_1_2_7
51 1_1_1_2_8 26 Cserum 0.000 3.000 1_1_1_2_8 1_1_1_2_8
52 1_1_1_2_8 26 Cserum 5.802 1.000 1_1_1_2_8 1_1_1_2_8
53 1_1_1_2_9 27 Cserum 0.000 3.100 1_1_1_2_9 1_1_1_2_9
54 1_1_1_2_9 27 Cserum 5.802 1.300 1_1_1_2_9 1_1_1_2_9
55 1_1_1_2_10 28 Cserum 0.000 1.500 1_1_1_2_10 1_1_1_2_10
56 1_1_1_2_10 28 Cserum 5.802 1.200 1_1_1_2_10 1_1_1_2_10
57 1_1_1_2_11 29 Cserum 0.000 3.100 1_1_1_2_11 1_1_1_2_11
58 1_1_1_2_11 29 Cserum 5.802 1.900 1_1_1_2_11 1_1_1_2_11
59 1_1_1_2_12 30 Cserum 0.000 1.900 1_1_1_2_12 1_1_1_2_12
60 1_1_1_2_12 30 Cserum 5.802 0.800 1_1_1_2_12 1_1_1_2_12
61 1_1_1_2_13 31 Cserum 0.000 1.000 1_1_1_2_13 1_1_1_2_13
62 1_1_1_2_13 31 Cserum 5.802 1.100 1_1_1_2_13 1_1_1_2_13
63 1_1_1_2_14 32 Cserum 0.000 1.200 1_1_1_2_14 1_1_1_2_14
64 1_1_1_2_14 32 Cserum 5.802 0.400 1_1_1_2_14 1_1_1_2_14
65 1_1_1_2_15 33 Cserum 0.000 2.500 1_1_1_2_15 1_1_1_2_15
66 1_1_1_2_15 33 Cserum 5.802 0.900 1_1_1_2_15 1_1_1_2_15
67 1_1_1_2_16 34 Cserum 0.000 1.200 1_1_1_2_16 1_1_1_2_16
68 1_1_1_2_16 34 Cserum 5.802 0.400 1_1_1_2_16 1_1_1_2_16
69 1_1_1_2_17 35 Cserum 0.000 2.300 1_1_1_2_17 1_1_1_2_17
70 1_1_1_2_17 35 Cserum 5.802 0.700 1_1_1_2_17 1_1_1_2_17
71 1_1_1_2_18 36 Cserum 0.000 1.500 1_1_1_2_18 1_1_1_2_18
72 1_1_1_2_18 36 Cserum 5.802 0.500 1_1_1_2_18 1_1_1_2_18
73 1_1_1_2_19 37 Cserum 0.000 1.900 1_1_1_2_19 1_1_1_2_19
74 1_1_1_2_19 37 Cserum 5.802 1.100 1_1_1_2_19 1_1_1_2_19
75 1_2_1 38 M_Cbgd_Css 2.200 5.710 1_2_1 1_2_1
76 1_3_1 39 M_Cbgd_Css 2.000 0.925 1_3_1 1_3_1
77 1_4_1 40 M_Cbgd_Css 2.000 1.060 1_4_1 1_4_1
dataset Sex City Train_Test datatype
1 Decatur M Train Male Decatur Train Individual
2 Decatur M Train Male Decatur Train Individual
3 Decatur M Train Male Decatur Train Individual
4 Decatur M Train Male Decatur Train Individual
5 Decatur M Train Male Decatur Train Individual
6 Decatur M Train Male Decatur Train Individual
7 Decatur M Train Male Decatur Train Individual
8 Decatur M Train Male Decatur Train Individual
9 Decatur M Train Male Decatur Train Individual
10 Decatur M Train Male Decatur Train Individual
11 Decatur M Train Male Decatur Train Individual
12 Decatur M Train Male Decatur Train Individual
13 Decatur M Train Male Decatur Train Individual
14 Decatur M Train Male Decatur Train Individual
15 Decatur M Train Male Decatur Train Individual
16 Decatur M Train Male Decatur Train Individual
17 Decatur M Train Male Decatur Train Individual
18 Decatur M Train Male Decatur Train Individual
19 Decatur F Train Female Decatur Train Individual
20 Decatur F Train Female Decatur Train Individual
21 Decatur F Train Female Decatur Train Individual
22 Decatur F Train Female Decatur Train Individual
23 Decatur F Train Female Decatur Train Individual
24 Decatur F Train Female Decatur Train Individual
25 Decatur F Train Female Decatur Train Individual
26 Decatur F Train Female Decatur Train Individual
27 Decatur F Train Female Decatur Train Individual
28 Decatur F Train Female Decatur Train Individual
29 Decatur F Train Female Decatur Train Individual
30 Decatur F Train Female Decatur Train Individual
31 Decatur F Train Female Decatur Train Individual
32 Decatur F Train Female Decatur Train Individual
33 Decatur F Train Female Decatur Train Individual
34 Decatur F Train Female Decatur Train Individual
35 Decatur F Train Female Decatur Train Individual
36 Decatur F Train Female Decatur Train Individual
37 Decatur M Test Male Decatur Test Individual
38 Decatur M Test Male Decatur Test Individual
39 Decatur M Test Male Decatur Test Individual
40 Decatur M Test Male Decatur Test Individual
41 Decatur M Test Male Decatur Test Individual
42 Decatur M Test Male Decatur Test Individual
43 Decatur M Test Male Decatur Test Individual
44 Decatur M Test Male Decatur Test Individual
45 Decatur M Test Male Decatur Test Individual
46 Decatur M Test Male Decatur Test Individual
47 Decatur M Test Male Decatur Test Individual
48 Decatur M Test Male Decatur Test Individual
49 Decatur M Test Male Decatur Test Individual
50 Decatur M Test Male Decatur Test Individual
51 Decatur M Test Male Decatur Test Individual
52 Decatur M Test Male Decatur Test Individual
53 Decatur M Test Male Decatur Test Individual
54 Decatur M Test Male Decatur Test Individual
55 Decatur F Test Female Decatur Test Individual
56 Decatur F Test Female Decatur Test Individual
57 Decatur F Test Female Decatur Test Individual
58 Decatur F Test Female Decatur Test Individual
59 Decatur F Test Female Decatur Test Individual
60 Decatur F Test Female Decatur Test Individual
61 Decatur F Test Female Decatur Test Individual
62 Decatur F Test Female Decatur Test Individual
63 Decatur F Test Female Decatur Test Individual
64 Decatur F Test Female Decatur Test Individual
65 Decatur F Test Female Decatur Test Individual
66 Decatur F Test Female Decatur Test Individual
67 Decatur F Test Female Decatur Test Individual
68 Decatur F Test Female Decatur Test Individual
69 Decatur F Test Female Decatur Test Individual
70 Decatur F Test Female Decatur Test Individual
71 Decatur F Test Female Decatur Test Individual
72 Decatur F Test Female Decatur Test Individual
73 Decatur F Test Female Decatur Test Individual
74 Decatur F Test Female Decatur Test Individual
75 Paulsboro-Train Mixed (all sexes) Paulsboro Train Summary
76 Horsham-Train Mixed (all sexes) Horsham Train Summary
77 Warminster-Test Mixed (all sexes) Warminster Test Summary
variable Dataset Time.desc Dataset.Time
1 Decatur M Train 1 Decatur M Train 1 T=0 Decatur M Train 1.T=0
2 Decatur M Train 1 Decatur M Train 1 T=5.802 Decatur M Train 1.T=5.802
3 Decatur M Train 2 Decatur M Train 2 T=0 Decatur M Train 2.T=0
4 Decatur M Train 2 Decatur M Train 2 T=5.802 Decatur M Train 2.T=5.802
5 Decatur M Train 3 Decatur M Train 3 T=0 Decatur M Train 3.T=0
6 Decatur M Train 3 Decatur M Train 3 T=5.802 Decatur M Train 3.T=5.802
7 Decatur M Train 4 Decatur M Train 4 T=0 Decatur M Train 4.T=0
8 Decatur M Train 4 Decatur M Train 4 T=5.802 Decatur M Train 4.T=5.802
9 Decatur M Train 5 Decatur M Train 5 T=0 Decatur M Train 5.T=0
10 Decatur M Train 5 Decatur M Train 5 T=5.802 Decatur M Train 5.T=5.802
11 Decatur M Train 6 Decatur M Train 6 T=0 Decatur M Train 6.T=0
12 Decatur M Train 6 Decatur M Train 6 T=5.802 Decatur M Train 6.T=5.802
13 Decatur M Train 7 Decatur M Train 7 T=0 Decatur M Train 7.T=0
14 Decatur M Train 7 Decatur M Train 7 T=5.802 Decatur M Train 7.T=5.802
15 Decatur M Train 8 Decatur M Train 8 T=0 Decatur M Train 8.T=0
16 Decatur M Train 8 Decatur M Train 8 T=5.802 Decatur M Train 8.T=5.802
17 Decatur M Train 9 Decatur M Train 9 T=0 Decatur M Train 9.T=0
18 Decatur M Train 9 Decatur M Train 9 T=5.802 Decatur M Train 9.T=5.802
19 Decatur F Train 10 Decatur F Train 10 T=0 Decatur F Train 10.T=0
20 Decatur F Train 10 Decatur F Train 10 T=5.802 Decatur F Train 10.T=5.802
21 Decatur F Train 11 Decatur F Train 11 T=0 Decatur F Train 11.T=0
22 Decatur F Train 11 Decatur F Train 11 T=5.802 Decatur F Train 11.T=5.802
23 Decatur F Train 12 Decatur F Train 12 T=0 Decatur F Train 12.T=0
24 Decatur F Train 12 Decatur F Train 12 T=5.802 Decatur F Train 12.T=5.802
25 Decatur F Train 13 Decatur F Train 13 T=0 Decatur F Train 13.T=0
26 Decatur F Train 13 Decatur F Train 13 T=5.802 Decatur F Train 13.T=5.802
27 Decatur F Train 14 Decatur F Train 14 T=0 Decatur F Train 14.T=0
28 Decatur F Train 14 Decatur F Train 14 T=5.802 Decatur F Train 14.T=5.802
29 Decatur F Train 15 Decatur F Train 15 T=0 Decatur F Train 15.T=0
30 Decatur F Train 15 Decatur F Train 15 T=5.802 Decatur F Train 15.T=5.802
31 Decatur F Train 16 Decatur F Train 16 T=0 Decatur F Train 16.T=0
32 Decatur F Train 16 Decatur F Train 16 T=5.802 Decatur F Train 16.T=5.802
33 Decatur F Train 17 Decatur F Train 17 T=0 Decatur F Train 17.T=0
34 Decatur F Train 17 Decatur F Train 17 T=5.802 Decatur F Train 17.T=5.802
35 Decatur F Train 18 Decatur F Train 18 T=0 Decatur F Train 18.T=0
36 Decatur F Train 18 Decatur F Train 18 T=5.802 Decatur F Train 18.T=5.802
37 Decatur M Test 19 Decatur M Test 19 T=0 Decatur M Test 19.T=0
38 Decatur M Test 19 Decatur M Test 19 T=5.802 Decatur M Test 19.T=5.802
39 Decatur M Test 20 Decatur M Test 20 T=0 Decatur M Test 20.T=0
40 Decatur M Test 20 Decatur M Test 20 T=5.802 Decatur M Test 20.T=5.802
41 Decatur M Test 21 Decatur M Test 21 T=0 Decatur M Test 21.T=0
42 Decatur M Test 21 Decatur M Test 21 T=5.802 Decatur M Test 21.T=5.802
43 Decatur M Test 22 Decatur M Test 22 T=0 Decatur M Test 22.T=0
44 Decatur M Test 22 Decatur M Test 22 T=5.802 Decatur M Test 22.T=5.802
45 Decatur M Test 23 Decatur M Test 23 T=0 Decatur M Test 23.T=0
46 Decatur M Test 23 Decatur M Test 23 T=5.802 Decatur M Test 23.T=5.802
47 Decatur M Test 24 Decatur M Test 24 T=0 Decatur M Test 24.T=0
48 Decatur M Test 24 Decatur M Test 24 T=5.802 Decatur M Test 24.T=5.802
49 Decatur M Test 25 Decatur M Test 25 T=0 Decatur M Test 25.T=0
50 Decatur M Test 25 Decatur M Test 25 T=5.802 Decatur M Test 25.T=5.802
51 Decatur M Test 26 Decatur M Test 26 T=0 Decatur M Test 26.T=0
52 Decatur M Test 26 Decatur M Test 26 T=5.802 Decatur M Test 26.T=5.802
53 Decatur M Test 27 Decatur M Test 27 T=0 Decatur M Test 27.T=0
54 Decatur M Test 27 Decatur M Test 27 T=5.802 Decatur M Test 27.T=5.802
55 Decatur F Test 28 Decatur F Test 28 T=0 Decatur F Test 28.T=0
56 Decatur F Test 28 Decatur F Test 28 T=5.802 Decatur F Test 28.T=5.802
57 Decatur F Test 29 Decatur F Test 29 T=0 Decatur F Test 29.T=0
58 Decatur F Test 29 Decatur F Test 29 T=5.802 Decatur F Test 29.T=5.802
59 Decatur F Test 30 Decatur F Test 30 T=0 Decatur F Test 30.T=0
60 Decatur F Test 30 Decatur F Test 30 T=5.802 Decatur F Test 30.T=5.802
61 Decatur F Test 31 Decatur F Test 31 T=0 Decatur F Test 31.T=0
62 Decatur F Test 31 Decatur F Test 31 T=5.802 Decatur F Test 31.T=5.802
63 Decatur F Test 32 Decatur F Test 32 T=0 Decatur F Test 32.T=0
64 Decatur F Test 32 Decatur F Test 32 T=5.802 Decatur F Test 32.T=5.802
65 Decatur F Test 33 Decatur F Test 33 T=0 Decatur F Test 33.T=0
66 Decatur F Test 33 Decatur F Test 33 T=5.802 Decatur F Test 33.T=5.802
67 Decatur F Test 34 Decatur F Test 34 T=0 Decatur F Test 34.T=0
68 Decatur F Test 34 Decatur F Test 34 T=5.802 Decatur F Test 34.T=5.802
69 Decatur F Test 35 Decatur F Test 35 T=0 Decatur F Test 35.T=0
70 Decatur F Test 35 Decatur F Test 35 T=5.802 Decatur F Test 35.T=5.802
71 Decatur F Test 36 Decatur F Test 36 T=0 Decatur F Test 36.T=0
72 Decatur F Test 36 Decatur F Test 36 T=5.802 Decatur F Test 36.T=5.802
73 Decatur F Test 37 Decatur F Test 37 T=0 Decatur F Test 37.T=0
74 Decatur F Test 37 Decatur F Test 37 T=5.802 Decatur F Test 37.T=5.802
75 Paulsboro-Train 38 Paulsboro-Train 38 T=2.2 Paulsboro-Train 38.T=2.2
76 Horsham-Train 39 Horsham-Train 39 T=2 Horsham-Train 39.T=2
77 Warminster-Test 40 Warminster-Test 40 T=2 Warminster-Test 40.T=2
#Multicheck plot
# Split Steady State Group into different populations for boxplot grouping
#df_check[df_check$Time.desc == "SteadyState" & grepl("Lubeck",df_check$Dataset),]$Time.desc <- "Lubeck"
#df_check[df_check$Time.desc == "SteadyState" & grepl("Little Hocking",df_check$Dataset),]$Time.desc <- "Little Hocking"
Modify aesthetics lookup table for boxplots
## additional source aesthetic lookup table for grey-scale time (years); merged legends save space on plotting output
times <- df_check%>% select(Time.desc, Time) %>% unique () %>%
mutate(rank = rank(Time) , grey = grey.colors(start=1,end=0.4, n = n()),
alpha = (rank)/8) %>%
select(-Time)
df_check <- df_check %>% mutate (legend_label = (paste0(City, "\n", Time.desc ) )) # add legend-labels
aes_lut <- df_check %>%
select(City, Train_Test, datatype,Time, Time.desc, legend_label) %>% unique () %>%
left_join(aes_lut[, c("City", "cols")], by = "City") %>% ungroup () %>% unique ()%>%
left_join (times, by = "Time.desc") %>%
arrange(datatype, City, Train_Test, Time) %>%
mutate(alpha = if_else(City == "Horsham", alpha/2, alpha)) %>% # otherwise too dark with this color
mutate_if(is.factor, as.character)
Changed grey start to 1 instead of 0.8, end at 0.6 instead of 0.4. Changed shape of symbols so they are filled.
##EF
df_decat <- df_check %>%
filter(City == "Decatur" & Train_Test %in% c ("Train", "Test")) %>%
mutate(panel = ordered (Train_Test, levels = c ("Train", "Test"),
labels = c("E: PFNA Decatur Train", "F: PFNA Decatur Test") ))
aes_lut_df_df_decat <- aes_lut %>%
filter(City == "Decatur" & Train_Test %in% c ("Train", "Test")) %>%
mutate_if(is.factor, as.character)
source( paste0(gsub(basename(here()), 'shared_functions', here()), '/plot_sum_boxplot.r'))
plt_train <- plot_sum_boxplot (dframe = df_decat, aes_lut= aes_lut_df_df_decat, facets = TRUE , pfas_nom = pfas_name )
print(plt_train)
ggsave(here ("output-plots",paste0( sa,"DecaturTrainTestboxplot",pfas_name,".pdf")),plt_train,dpi=600)
Saving 6.5 x 3.5 in image
ggsave(here ("output-plots",paste0( sa,"DecaturTrainTestboxplot",pfas_name,".png")),plt_train,dpi=600)
Saving 6.5 x 3.5 in image
Changed grey start to 1 instead of 0.8, end at 0.6 instead of 0.4. Added shapes and fills to data points.
lets <- LETTERS;
names(lets)[1:(length(unique(df_check$dataset))-4)]<-as.character(unique(df_check$dataset))[5:length(unique(df_check$dataset))]
for (d in unique(df_check$dataset)) { # d = unique(df_check$dataset)[11]
ddset <- df_check %>%
filter(dataset == d)
aes_lut_ddset <- ddset %>% select(legend_label, City,Train_Test,datatype, Time.desc ) %>% unique () %>% inner_join(aes_lut)
gt <- ifelse(is.na(lets[d]),d,paste0(lets[d],": ", d))
plt <- plot_sum_boxplot(dframe = ddset, aes_lut= aes_lut_ddset, gtitle= gt, facets = FALSE, pfas_nom = pfas_name)
print(plt)
ggsave(here ("output-plots",
paste0( sa, d,"-boxplot-",
pfas_name,".pdf")) ,
plt,dpi=600)
ggsave(here ("output-plots",
paste0( sa, d,"-boxplot-",
pfas_name,".png")) ,
plt,dpi=600)
}
Joining, by = c("legend_label", "City", "Train_Test", "datatype", "Time.desc")
Saving 6.5 x 3.5 in image
Saving 6.5 x 3.5 in image
Joining, by = c("legend_label", "City", "Train_Test", "datatype", "Time.desc")
Saving 6.5 x 3.5 in image
Saving 6.5 x 3.5 in image
Joining, by = c("legend_label", "City", "Train_Test", "datatype", "Time.desc")
Saving 6.5 x 3.5 in image
Saving 6.5 x 3.5 in image
Joining, by = c("legend_label", "City", "Train_Test", "datatype", "Time.desc")
Saving 6.5 x 3.5 in image
Saving 6.5 x 3.5 in image
Joining, by = c("legend_label", "City", "Train_Test", "datatype", "Time.desc")
Saving 6.5 x 3.5 in image
Saving 6.5 x 3.5 in image
Joining, by = c("legend_label", "City", "Train_Test", "datatype", "Time.desc")
Saving 6.5 x 3.5 in image
Saving 6.5 x 3.5 in image
Joining, by = c("legend_label", "City", "Train_Test", "datatype", "Time.desc")
Saving 6.5 x 3.5 in image
Saving 6.5 x 3.5 in image
### make Training plot
df_d_trt <- df_check %>%
filter( (Train_Test == "Train") & ((Output_Var == "M_Cbgd_Css") | (Output_Var == "M_Cserum"))) %>%
mutate_if(is.factor, as.character) %>% # drop factor levels unused
mutate(Dataset.Time = factor(Dataset.Time))
aes_lut_df_d_trt <- df_d_trt %>% select(City, datatype,Time, Time.desc, legend_label) %>%
inner_join(aes_lut ) %>%
select(-Train_Test) %>% ungroup () %>% unique ()
Joining, by = c("City", "datatype", "Time", "Time.desc", "legend_label")
plt_train <- plot_sum_boxplot(dframe = df_d_trt, aes_lut= aes_lut_df_d_trt,
gtitle="A: Summary Data - Train" , facets = FALSE, pfas_nom = pfas_name )
print(plt_train)
ggsave(here ("output-plots", paste0( sa, "SummaryTrainDataboxplot",pfas_name,".pdf")), plt_train,dpi=600)
Saving 6.5 x 3.5 in image
ggsave(here ("output-plots", paste0( sa, "SummaryTrainDataboxplot",pfas_name,".png")), plt_train,dpi=600)
Saving 6.5 x 3.5 in image
### make Test plot
df_d_test <- df_check %>%
filter((Train_Test == "Test") &
((Output_Var == "M_Cbgd_Css") | (Output_Var == "M_Cserum"))) %>%
mutate_if(is.factor, as.character) %>% # drop factor levels unused
mutate(Dataset.Time = factor(Dataset.Time))
aes_lut_df_d_test <- df_d_test %>% select(City, datatype,Time, Time.desc, legend_label) %>%
inner_join(aes_lut ) %>%
select(-Train_Test) %>% ungroup () %>% unique ()
Joining, by = c("City", "datatype", "Time", "Time.desc", "legend_label")
plt_test <- plot_sum_boxplot(dframe = df_d_test, aes_lut= aes_lut_df_d_test,
gtitle="B: Summary Data - Test", facets = FALSE, pfas_nom = pfas_name)
print(plt_test)
ggsave(here ("output-plots",paste0( sa, "SummaryTestDataboxplot",pfas_name,".pdf")), plt_test,dpi=600)
Saving 6.5 x 3.5 in image
ggsave(here ("output-plots",paste0( sa, "SummaryTestDataboxplot",pfas_name,".png")), plt_test,dpi=600)
Saving 6.5 x 3.5 in image
Shows shift in background estimate.
gmscale<-0.8
dat <- multicheck$parms.samp[,grep("M_ln_Cbgd",names(multicheck$parms.samp))]
datasetnames <- as.character(unique(calibdata$dataset))
datasetnames <- gsub(" Train","-Train",datasetnames)
datasetnames <- gsub(" Test","-Train",datasetnames)
datasetnames <- gsub(" M","",datasetnames)
datasetnames <- gsub(" F","",datasetnames)
datasetnames<-datasetnames[!duplicated(datasetnames)]
names(dat) <- datasetnames
dat <- dat[,grep("Train",names(dat))]
dat.df <- pivot_longer(dat,1:ncol(dat))
dat.df <- rbind(dat.df,
data.frame(name="Prior",value=rnorm(5000,m=log(gmscale),sd=0.4055)))
dat.df$name <- factor(dat.df$name,levels=rev(
c("Prior",datasetnames[grep("Train",datasetnames)])))
dat.df$value <- exp(dat.df$value)
p<-ggplot(dat.df)+
#geom_violin(aes(x=name,y=value,fill=name=="Prior"))+
geom_boxplot(aes(x=name,y=value,fill=name=="Prior"),outlier.shape=NA)+
scale_y_log10()+coord_flip()+
scale_fill_manual(name=NULL,
values=c("#009988", "#EE7733" ))+
theme_classic() +
geom_hline(yintercept = gmscale,color="grey")+
theme(legend.position="none",
panel.background = element_rect(color="black",size=1))+
ylab("Posterior shift in Background Concentration")
print(p)
ggsave(here ("output-plots",paste0( sa, "PFNA_GM_Cbgd.pdf")) ,p,dpi=600)
Saving 5 x 6 in image
ggsave(here ("output-plots",paste0( sa, "PFNA_GM_Cbgd.png")) ,p,dpi=600)
Saving 5 x 6 in image
For PFNA, the population GM of the half-life has a posterior distribution that is narrower than the prior, with a posterior median (95% CI) estimate of 3.06 (2.16-4.37) years. The population GSD posterior is larger than the prior at 1.47(1.44-1.75).
dat <- multicheck$parms.samp[,c("M_ln_k.1.","V_ln_k.1.", "M_ln_Vd.1.", "SD_ln_Vd.1.")]
names(dat) <- c("M_ln_k(1)","V_ln_k(1)", "M_ln_Vd(1)", "SD_ln_Vd(1)")
set.seed(3.14159)
dat$z_ln_k <- rnorm(nrow(dat))
dat$z_ln_Vd <- rnorm(nrow(dat))
dat %>% rename_()
dat$ln_k_i <- dat$`M_ln_k(1)` + sqrt(dat$`V_ln_k(1)`)*dat$z_ln_k
dat$ln_Vd_i <- dat$`M_ln_Vd(1)`+ dat$`SD_ln_Vd(1)`*dat$z_ln_Vd
linmod <- lm(ln_Vd_i ~ ln_k_i,data=dat)
ggplot(dat) + geom_point(aes(ln_k_i,ln_Vd_i)) +
labs(subtitle=paste("Adj R2 =",signif(summary(linmod)$adj.r.squared,2)))
qqnorm(dat$ln_k_i,main="ln k Q-Q Normal")
qqline(dat$ln_k_i,col="red")
plot(ecdf(dat$ln_k_i))
x <- seq(-3,1,0.01)
m_ln_k_i <- mean(dat$ln_k_i)
sd_ln_k_i <- sd(dat$ln_k_i)
lines(x,pnorm(x,mean=m_ln_k_i,sd=sd_ln_k_i),col="red")
text(m_ln_k_i-2*sd_ln_k_i,0.9,paste("m =",signif(m_ln_k_i,4),"\nsd =",signif(sd_ln_k_i,4)))
qqnorm(dat$ln_Vd_i,main="ln Vd Q-Q Normal")
qqline(dat$ln_Vd_i,col="red")
plot(ecdf(dat$ln_Vd_i))
x <- seq(-3,1,0.01)
m_ln_Vd_i <- mean(dat$ln_Vd_i)
sd_ln_Vd_i <- sd(dat$ln_Vd_i)
lines(x,pnorm(x,mean=m_ln_Vd_i,sd=sd_ln_Vd_i),col="red")
text(m_ln_Vd_i-2*sd_ln_Vd_i,0.9,paste("m =",signif(m_ln_Vd_i,4),"\nsd =",signif(sd_ln_Vd_i,4)))
hl_i <- log(2)/ exp(dat$ln_k_i) # individual half-life
med_hl_i <- paste(signif (median (hl_i), 3)) # median of individual half-life
ci_med_hl_i <- paste(signif (quantile(hl_i, prob=c(0.025,0.975)), 3),collapse="-") # 95ci med individual halflife
ci98_med_hl_i <- paste(signif (quantile(hl_i, prob=c(0.01,0.99)), 3),collapse="-") # 98ci med individual halflife
gm_hl_i <- paste(signif (exp(mean(log(hl_i))), 3)) # gm (which should be really close)
gsd_hl_i <- paste(signif (exp(sd(log(hl_i))), 3)) # gsd individual
med_Vd_i <- paste(signif (median(exp(dat$ln_Vd_i)), 3)) # median individual Vd
ci_med_Vd_i <-paste(signif (quantile(exp(dat$ln_Vd_i), prob=c(0.025,0.975)), 3),collapse="-") # 95ci med individual Vd
ci98_med_Vd_i <-paste(signif (quantile(exp(dat$ln_Vd_i), prob=c(0.01,0.99)), 3),collapse="-") # 98ci med individual Vd
gm_vd_i <- paste(signif (exp(mean(dat$ln_Vd_i)), 3)) # gm (which should be really close)
gsd_vd_i<- paste(signif (exp(sd(dat$ln_Vd_i)), 3)) # gsd indiv
med_CL_i <- paste(signif (median(exp(dat$ln_Vd_i+dat$ln_k_i)), 3)) # median individual CL
ci_med_CL_i <-paste(signif (quantile(exp(dat$ln_Vd_i+dat$ln_k_i), prob=c(0.025,0.975)), 3),collapse="-") # 95ci med individual CL
ci98_med_CL_i <-paste(signif (quantile(exp(dat$ln_Vd_i+dat$ln_k_i), prob=c(0.01,0.99)), 3),collapse="-") # 98ci med individual CL
gm_CL_i <- paste(signif (exp(mean(dat$ln_Vd_i+dat$ln_k_i)), 3)) # gm (which should be really close)
gsd_CL_i<- paste(signif (exp(sd(dat$ln_Vd_i+dat$ln_k_i)), 3)) # gsd indiv
PFNA_priors <- data.frame(
halflife_GM= log(2)/rlnorm(50000,
meanlog=-1.80181,sdlog=0.4055))
M_k <- exp(as.numeric(dat$`M_ln_k(1)`))
PFNA_halflife_GM <- log(2)/M_k
PFNA_hlgm_pr_med <- signif(median(PFNA_priors$halflife_GM,3))
PFNA_hlgm_pr_med_95ci <-paste(signif(quantile(PFNA_priors$halflife_GM,
prob=c(0.025,0.975)),
3),
collapse="-")
PFNA_hl_median_gm <- signif(median(PFNA_halflife_GM),3)
PFNA_hl_median_gm_95ci <- paste(signif(quantile(PFNA_halflife_GM,
prob=c(0.025,0.975)),3),collapse="-")
p<-ggplot()+
stat_density(aes(halflife_GM, color = "Prior"),data=PFNA_priors,geom="line",size=2)+
stat_density(aes(PFNA_halflife_GM,stat(density),color="Posterior"),geom="line",size=1.5)+
xlim(0,15)+
labs(title = bquote("E: PFNA"~T[1/2]~"Population GM") ,
subtitle=paste("Posterior Median (95% CI): ",
PFNA_hl_median_gm," (",
PFNA_hl_median_gm_95ci,
")",sep=""))+
xlab(bquote("Population GM"~T[1/2]~"(yrs)")) +
scale_color_manual(name=NULL,
values=c(Prior="#009988", Posterior="#EE7733" )) +
theme_classic() +
theme(legend.title = element_blank(),legend.position=c(0.8,0.7),
panel.background = element_rect(color="black",size=1),
legend.background = element_rect(fill="transparent", color=NA))
print(p)
Warning: Removed 41 rows containing non-finite values (stat_density).
ggsave(here ("output-plots", paste0( sa, "PFNA_hl_gm.pdf")) ,p,dpi=600)
Saving 4 x 2.5 in image
Warning: Removed 41 rows containing non-finite values (stat_density).
ggsave(here ("output-plots", paste0( sa, "PFNA_hl_gm.png")) ,p,dpi=600)
Saving 4 x 2.5 in image
Warning: Removed 41 rows containing non-finite values (stat_density).
PFNA_priors$halflife_GSD = exp(sqrt(exp(rnorm(50000,m=log(0.2000),sd=log(1.275)))))
PFNA_halflife_GSD <- exp(sqrt(dat$`V_ln_k(1)`))
PFNA_hlgsd_pr_med <- signif(median(PFNA_priors$halflife_GSD,3))
PFNA_hlgsd_pr_med_95ci <-paste(signif(quantile(PFNA_priors$halflife_GSD,
prob=c(0.025,0.975)),
3),
collapse="-")
PFNA_hl_gsd_med <- signif(median(PFNA_halflife_GSD),3)
PFNA_hl_gsd_med_95ci <- paste(signif(quantile(PFNA_halflife_GSD,
prob=c(0.025,0.975)),3),collapse="-")
p<-ggplot()+
stat_density(aes(halflife_GSD, color = "Prior"),data=PFNA_priors,geom="line",size=2)+
stat_density(aes(PFNA_halflife_GSD,stat(density), color = "Posterior"),geom="line",size=1.5)+
xlim(1,3)+
labs(title = bquote("F: PFNA"~T[1/2]~"Population GSD"),
subtitle=paste("Posterior Median (95% CI): ",
PFNA_hl_gsd_med," (",
PFNA_hl_gsd_med_95ci,
")",sep=""))+
xlab(bquote("Population GSD"~T[1/2]))+
scale_color_manual(name=NULL,
values=c(Prior="#009988", Posterior="#EE7733" )) +
theme_classic() +
theme(legend.title = element_blank(),legend.position=c(0.8,0.7),
panel.background = element_rect(color="black",size=1),
legend.background = element_rect(fill="transparent", color=NA))
print(p)
ggsave(here ("output-plots", paste0( sa, "PFNA_hl_gsd.pdf")), p,dpi=600)
ggsave(here ("output-plots", paste0( sa, "PFNA_hl_gsd.png")), p,dpi=600)
For PFNA, the data were not particularly informative, but slightly increased the estimate of the median to 0.308(0.223-0.548) slightly. They were not informative as to the population GSD, with the posterior distributions essentially unchanged from the priors.
PFNA_priors$Vd_GM <- rlnorm(50000,
meanlog=-1.77196,
sdlog=0.2624)
PFNA_Vd_GM <- exp(dat$`M_ln_Vd(1)`)
PFNA_vd_gm_pr_med <- signif(median(PFNA_priors$Vd_GM,3))
PFNA_vd_gm_pr_med_95ci <- paste(signif(quantile(PFNA_priors$Vd_GM,
prob=c(0.025,0.975)), 3), collapse="-")
PFNA_vd_gm_med <- signif(median(PFNA_Vd_GM),3)
PFNA_vd_gm_med_95ci <- paste(signif(quantile(PFNA_Vd_GM,
prob=c(0.025,0.975)),3),collapse="-")
p<-ggplot()+
stat_density(aes(Vd_GM, color = "Prior"),data=PFNA_priors,geom="line",size=2)+
stat_density(aes(PFNA_Vd_GM,stat(density), color = "Posterior"),geom="line",size=1.5)+
xlim(0,1)+labs(title = bquote("E: PFNA"~V[d]~"Population GM"),
subtitle=paste("Posterior Median (95% CI): ",
PFNA_vd_gm_med," (",
PFNA_vd_gm_med_95ci,")",sep=""))+
xlab(bquote("Population GM"~V[d]~"(l/kg)"))+
scale_color_manual(name=NULL,
values=c(Prior="#009988", Posterior="#EE7733" )) +
theme_classic() +
theme(legend.title = element_blank(),legend.position=c(0.8,0.7),
panel.background = element_rect(color="black",size=1),
legend.background = element_rect(fill="transparent", color=NA))
print(p)
ggsave(here ("output-plots",paste0( sa, "PFNA_vd_gm.pdf")), p,dpi=600)
ggsave(here ("output-plots",paste0( sa, "PFNA_vd_gm.png")), p,dpi=600)
PFNA_priors$Vd_GSD = exp(abs(rnorm(50000,sd=0.17)))
PFNA_Vd_GSD <- exp(dat$`SD_ln_Vd(1)`)
PFNA_vd_gsd_pr_med <- signif(median(PFNA_priors$Vd_GSD,3))
PFNA_vd_gsd_pr_med_95ci <- paste(signif(quantile(PFNA_priors$Vd_GSD,
prob=c(0.025,0.975)), 3), collapse="-")
PFNA_vd_gsd_med <- signif(median(PFNA_Vd_GSD),3)
PFNA_vd_gsd_med_95ci <- paste(signif(quantile(PFNA_Vd_GSD,
prob=c(0.025,0.975)),3),collapse="-")
p<-ggplot()+
stat_density(aes(Vd_GSD, color = "Prior"),data=PFNA_priors,geom="line",size=2)+
stat_density(aes(PFNA_Vd_GSD,stat(density), color = "Posterior"),geom="line",size=1.5)+
xlim(1,3)+
labs(title = bquote("F: PFNA"~V[d]~"Population GSD "),
subtitle=paste("Posterior Median (95% CI): ",
PFNA_vd_gsd_med," (",
PFNA_vd_gsd_med_95ci,
")",sep=""))+
xlab(bquote("Population GSD"~V[d]))+
scale_color_manual(name=NULL,
values=c(Prior="#009988", Posterior="#EE7733" )) +
theme_classic() +
theme(legend.title = element_blank(),legend.position=c(0.8,0.7),
panel.background = element_rect(color="black",size=1),
legend.background = element_rect(fill="transparent", color=NA))
print(p)
ggsave(here ("output-plots",paste0( sa, "PFNA_vd_gsd.pdf")), p,dpi=600)
ggsave(here ("output-plots",paste0( sa, "PFNA_vd_gsd.png")), p,dpi=600)
Cl is k * Vd
PFNA_priors$CL_GM <- PFNA_priors$Vd_GM * (log(2)/PFNA_priors$halflife_GM)
PFNA_CL_GM <- exp(dat$`M_ln_Vd(1)` + dat$`M_ln_k(1)`)
PFNA_cl_gm_pr_med <- signif(median(PFNA_priors$CL_GM,3))
PFNA_cl_gm_pr_med_95ci <- paste(signif(quantile(PFNA_priors$CL_GM,
prob=c(0.025,0.975)), 3), collapse="-")
PFNA_cl_gm_med <- signif(median(PFNA_CL_GM),3)
PFNA_cl_gm_med_95ci <- paste(signif(quantile(PFNA_CL_GM,
prob=c(0.025,0.975)),3),collapse="-")
p<-ggplot()+
stat_density(aes(CL_GM, color = "Prior"),data=PFNA_priors,geom="line",size=2)+
stat_density(aes(PFNA_CL_GM,stat(density), color = "Posterior"),geom="line",size=1.5)+
xlim(0,0.25)+labs(title = "C: PFNA Clearance Pop. GM ",subtitle=paste("Posterior Median (95% CI): ",
PFNA_cl_gm_med," (",
PFNA_cl_gm_med_95ci,
")",sep=""))+
xlab("Pop. GM CL (l/(kg-yr))")+
scale_color_manual(name=NULL,
values=c(Prior="#009988", Posterior="#EE7733" )) +
theme_classic() +
theme(legend.title = element_blank(),legend.position=c(0.8,0.7),
panel.background = element_rect(color="black",size=1),
legend.background = element_rect(fill="transparent", color=NA))
print(p)
ggsave(here ("output-plots",paste0( sa, "PFNA_CL_gm.pdf")) ,p,dpi=600)
ggsave(here ("output-plots",paste0( sa, "PFNA_CL_gm.png")) ,p,dpi=600)
PFNA_priors$CL_GSD = exp(sqrt(log(PFNA_priors$Vd_GSD)^2 +
log(PFNA_priors$halflife_GSD)^2))
PFNA_CL_GSD <- exp(sqrt(log(PFNA_Vd_GSD)^2 +
log(PFNA_halflife_GSD)^2))
PFNA_CL_gsd_pr_med <- signif(median(PFNA_priors$CL_GSD,3))
PFNA_CL_gsd_pr_med_95ci <- paste(signif(quantile(PFNA_priors$CL_GSD,
prob=c(0.025,0.975)), 3), collapse="-")
PFNA_CL_gsd_med <- signif(median(PFNA_CL_GSD),3)
PFNA_CL_gsd_med_95ci <- paste(signif(quantile(PFNA_CL_GSD,
prob=c(0.025,0.975)),3),collapse="-")
p<-ggplot()+
stat_density(aes(CL_GSD, color = "Prior"),data=PFNA_priors,geom="line",size=2)+
stat_density(aes(PFNA_CL_GSD,stat(density), color = "Posterior"),geom="line",size=1.5)+
xlim(1,3)+
labs(title = bquote("H: PFNA"~CL~"Population GSD "),
subtitle=paste("Posterior Median (95% CI): ",
PFNA_CL_gsd_med," (",
PFNA_CL_gsd_med_95ci,
")",sep=""))+
xlab(bquote("Population GSD"~CL))+
scale_color_manual(name=NULL,
values=c(Prior="#009988", Posterior="#EE7733" )) +
theme_classic() +
theme(legend.title = element_blank(),legend.position=c(0.8,0.7),
panel.background = element_rect(color="black",size=1),
legend.background = element_rect(fill="transparent", color=NA))
print(p)
ggsave(here ("output-plots",paste0( sa,"PFNA_CL_gsd.pdf")) ,p,dpi=600)
ggsave(here ("output-plots",paste0( sa,"PFNA_CL_gsd.png")) ,p,dpi=600)
PFNA_hlgm_pr_med <- paste(signif(PFNA_hlgm_pr_med, 3))
PFNA_hl_median_gm<- paste(signif(PFNA_hl_median_gm, 3))
PFNA_hlgsd_pr_med<- paste(signif(PFNA_hlgsd_pr_med, 3))
PFNA_hl_gsd_med<- paste(signif(PFNA_hl_gsd_med, 3))
PFNA_vd_gm_pr_med<- paste(signif(PFNA_vd_gm_pr_med, 3))
PFNA_vd_gm_med<- paste(signif(PFNA_vd_gm_med, 3))
PFNA_vd_gsd_pr_med<- paste(signif(PFNA_vd_gsd_pr_med, 3))
PFNA_vd_gsd_med<- paste(signif(PFNA_vd_gsd_med, 3))
PFNA_cl_gm_pr_med<- paste(signif(PFNA_cl_gm_pr_med, 3))
PFNA_cl_gm_med<- paste(signif(PFNA_cl_gm_med, 3))
| Parameter | Prior GM | Posterior GM | Prior GSD | Posterior GSD |
|---|---|---|---|---|
| Half-life (years) | 4.2 | 2.35 | 1.56 | 1.53 |
| HL [95% CI] | [1.89-9.37] | [1.65-3.16] | [1.42-1.77] | [1.4-1.7] |
| Volume of distribution | 0.17 | 0.186 | 1.12 | 1.12 |
| \(V_D\) [95% CI] | [0.102-0.284] | [0.113-0.302] | [1.01-1.46] | [1.01-1.51] |
| Clearance | 0.028 | 0.0559 | ||
| \(CL\) [95% CI] | [0.0108-0.0727] | [0.0327-0.0932] | [] | [] |
| Parameter | median GM [95% CI] [[98% CI]] | GM calculator input | GSD individual |
|---|---|---|---|
| Half-life (years) | 2.27 [ 0.826-5.36 ] [[ 0.758-5.94 ]] | 2.27 | 1.61 |
| Volume of distribution \(V_D\) | 0.183 [ 0.0998-0.324 ] [[ 0.085-0.398 ]] | 0.185 | 1.35 |
| Clearance (L/kg-yr) | 0.0573 [ 0.0188-0.163 ] [[ 0.0165-0.199 ]] | 0.0563 | 1.7 |
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 3.6.3 (2020-02-29)
os Red Hat Enterprise Linux Server 7.9 (Maipo)
system x86_64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/New_York
date 2022-01-23
─ Packages ───────────────────────────────────────────────────────────────────
package * version date lib source
assertthat 0.2.1 2019-03-21 [2] CRAN (R 3.6.3)
backports 1.2.1 2020-12-09 [2] CRAN (R 3.6.3)
bayesplot * 1.8.0 2021-01-10 [2] CRAN (R 3.6.3)
broom 0.7.5 2021-02-19 [2] CRAN (R 3.6.3)
bslib 0.2.4 2021-01-25 [2] CRAN (R 3.6.3)
cachem 1.0.4 2021-02-13 [2] CRAN (R 3.6.3)
callr 3.5.1 2020-10-13 [2] CRAN (R 3.6.3)
cellranger 1.1.0 2016-07-27 [2] CRAN (R 3.6.3)
cli 2.3.1 2021-02-23 [2] CRAN (R 3.6.3)
coda * 0.19-4 2020-09-30 [2] CRAN (R 3.6.3)
codetools 0.2-18 2020-11-04 [2] CRAN (R 3.6.3)
colorspace 2.0-0 2020-11-11 [2] CRAN (R 3.6.3)
crayon 1.4.1 2021-02-08 [2] CRAN (R 3.6.3)
DBI 1.1.1 2021-01-15 [2] CRAN (R 3.6.3)
dbplyr 2.1.0 2021-02-03 [2] CRAN (R 3.6.3)
debugme 1.1.0 2017-10-22 [2] CRAN (R 3.6.3)
desc 1.3.0 2021-03-05 [2] CRAN (R 3.6.3)
devtools 2.3.2 2020-09-18 [2] CRAN (R 3.6.3)
digest 0.6.27 2020-10-24 [2] CRAN (R 3.6.3)
dplyr * 1.0.5 2021-03-05 [2] CRAN (R 3.6.3)
ellipsis 0.3.1 2020-05-15 [2] CRAN (R 3.6.3)
evaluate 0.14 2019-05-28 [2] CRAN (R 3.6.3)
fansi 0.4.2 2021-01-15 [2] CRAN (R 3.6.3)
farver 2.1.0 2021-02-28 [2] CRAN (R 3.6.3)
fastmap 1.1.0 2021-01-25 [2] CRAN (R 3.6.3)
forcats * 0.5.1 2021-01-27 [2] CRAN (R 3.6.3)
fs 1.5.0 2020-07-31 [2] CRAN (R 3.6.3)
generics 0.1.0 2020-10-31 [2] CRAN (R 3.6.3)
ggplot2 * 3.3.3 2020-12-30 [2] CRAN (R 3.6.3)
ggridges 0.5.3 2021-01-08 [2] CRAN (R 3.6.3)
ggsci * 2.9 2018-05-14 [2] CRAN (R 3.6.3)
glue 1.4.2 2020-08-27 [2] CRAN (R 3.6.3)
gtable 0.3.0 2019-03-25 [2] CRAN (R 3.6.3)
haven 2.3.1 2020-06-01 [2] CRAN (R 3.6.3)
here * 1.0.1 2020-12-13 [2] CRAN (R 3.6.3)
highr 0.8 2019-03-20 [2] CRAN (R 3.6.3)
hms 1.0.0 2021-01-13 [2] CRAN (R 3.6.3)
htmltools 0.5.1.1 2021-01-22 [2] CRAN (R 3.6.3)
httr 1.4.2 2020-07-20 [2] CRAN (R 3.6.3)
jquerylib 0.1.3 2020-12-17 [2] CRAN (R 3.6.3)
jsonlite 1.7.2 2020-12-09 [2] CRAN (R 3.6.3)
khroma * 1.7.0 2021-09-02 [1] CRAN (R 3.6.3)
knitr 1.31 2021-01-27 [2] CRAN (R 3.6.3)
labeling 0.4.2 2020-10-20 [2] CRAN (R 3.6.3)
lattice 0.20-41 2020-04-02 [2] CRAN (R 3.6.3)
lifecycle 1.0.0 2021-02-15 [2] CRAN (R 3.6.3)
lubridate 1.7.10 2021-02-26 [2] CRAN (R 3.6.3)
magrittr 2.0.1 2020-11-17 [2] CRAN (R 3.6.3)
memoise 2.0.0 2021-01-26 [2] CRAN (R 3.6.3)
modelr 0.1.8 2020-05-19 [2] CRAN (R 3.6.3)
munsell 0.5.0 2018-06-12 [2] CRAN (R 3.6.3)
pillar 1.5.1 2021-03-05 [2] CRAN (R 3.6.3)
pkgbuild 1.2.0 2020-12-15 [2] CRAN (R 3.6.3)
pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 3.6.3)
pkgload 1.2.0 2021-02-23 [2] CRAN (R 3.6.3)
plyr 1.8.6 2020-03-03 [2] CRAN (R 3.6.3)
prettyunits 1.1.1 2020-01-24 [2] CRAN (R 3.6.3)
pROC 1.17.0.1 2021-01-13 [2] CRAN (R 3.6.3)
processx 3.4.5 2020-11-30 [2] CRAN (R 3.6.3)
ps 1.6.0 2021-02-28 [2] CRAN (R 3.6.3)
purrr * 0.3.4 2020-04-17 [2] CRAN (R 3.6.3)
R6 2.5.0 2020-10-28 [2] CRAN (R 3.6.3)
Rcpp 1.0.6 2021-01-15 [2] CRAN (R 3.6.3)
readr * 1.4.0 2020-10-05 [2] CRAN (R 3.6.3)
readxl 1.3.1 2019-03-13 [2] CRAN (R 3.6.3)
remotes 2.2.0 2020-07-21 [2] CRAN (R 3.6.3)
reprex 1.0.0 2021-01-27 [2] CRAN (R 3.6.3)
reshape2 * 1.4.4 2020-04-09 [2] CRAN (R 3.6.3)
rlang 0.4.10 2020-12-30 [2] CRAN (R 3.6.3)
rmarkdown 2.7 2021-02-19 [2] CRAN (R 3.6.3)
rprojroot 2.0.2 2020-11-15 [2] CRAN (R 3.6.3)
rstudioapi 0.13 2020-11-12 [2] CRAN (R 3.6.3)
rvest 1.0.0 2021-03-09 [2] CRAN (R 3.6.3)
sass 0.3.1 2021-01-24 [2] CRAN (R 3.6.3)
scales * 1.1.1 2020-05-11 [2] CRAN (R 3.6.3)
sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 3.6.3)
stringi 1.5.3 2020-09-09 [2] CRAN (R 3.6.3)
stringr * 1.4.0 2019-02-10 [2] CRAN (R 3.6.3)
testthat 3.0.2 2021-02-14 [2] CRAN (R 3.6.3)
tibble * 3.1.0 2021-02-25 [2] CRAN (R 3.6.3)
tidyr * 1.1.3 2021-03-03 [2] CRAN (R 3.6.3)
tidyselect 1.1.0 2020-05-11 [2] CRAN (R 3.6.3)
tidyverse * 1.3.0 2019-11-21 [2] CRAN (R 3.6.3)
usethis 2.0.1 2021-02-10 [2] CRAN (R 3.6.3)
utf8 1.2.1 2021-03-12 [2] CRAN (R 3.6.3)
vctrs 0.3.6 2020-12-17 [2] CRAN (R 3.6.3)
withr 2.4.1 2021-01-26 [2] CRAN (R 3.6.3)
xfun 0.22 2021-03-11 [2] CRAN (R 3.6.3)
xml2 1.3.2 2020-04-23 [2] CRAN (R 3.6.3)
yaml 2.2.1 2020-02-01 [2] CRAN (R 3.6.3)
yardstick * 0.0.9 2021-11-22 [1] CRAN (R 3.6.3)
[1] /home/ad.abt.local/layc/R/x86_64-pc-linux-gnu-library/3.6
[2] /opt/R/3.6.3/lib64/R/library